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ABSTRACT 

We argue that type la supernovae (SNe la) are the result of head-on collisions of White Dwarfs (WDs) in 
triple systems. The thermonuclear explosions resulting from the zero-impact-parameter collisions of WDs 
are calculated from first principles by using 2D hydrodynamical simulations. Collisions of typical WDs with 
masses 0.5-0.9M Q result in explosions that synthesize 56 Ni masses in the range of 0.15-0.8M Q , spanning the 
wide distribution of yields observed for the majority of SNe la. The robustness of the shock ignition process is 
verified with a detailed study using a one-dimensional toy model and analytic tools. The late-time (> 50 days 
after peak) bolometric light curve is equal to the instantaneous energy deposition and is calculated exactly, by 
solving the transport of 7-rays emitted by the decay of 56 Ni using a Monte-Carlo code. All collisions are found 
to have the same late-time light curves, when normalized to the amount of synthesized Ni. This universal 
light curve is shown to agree with the majority of the supernovae in the compilation made by M. Stritzinger to 
an accuracy of better than 30% in the range 40 < t < 80 days after bolometric peak. The widths of the 56 M- 
mass-weighted-line-of-sight velocity distributions are correlated with the 5b Ni yield and in agreement with the 
observed Mazzali relation. The continuous distribution of observed SN la features, is naturally reproduced with 
the distribution of WD masses involved in the collisions. The effect of a non-zero impact parameter requires 
further studies, using 3D hydrodynamical codes. 

Subject headings: hydrodynamics - shock waves - supernovae: individual (la) 



The explosion mechanism of carbon-oxygen white dwarfs 
(CO WDs) to produce type la SNe is unknown (see 
iHillebrandt & Niemeverl 120001 for a review). It is widely 
thought that the explosion is mediated by accretion of mat- 
ter onto the WD. As it approaches the Chandrasekhar mass, 
an unknown process is believed to ignite a nuclear detonation. 
The observed light curve is known to be powered by 56 M syn- 
thesized in the explosion through the radioactive decay chain 
56 Ni -> 56 Co ^ 6 Fe. 

iKatz & Dong! d2012h suggested that some or all of SNe la 
arise from head-on collisions of typical field WDs in triple 
systems. In this scenario, the nuclear detonation is due to 
a shock ignition and is not related to whether or not the to- 
tal mass of the colliding WDs is above or below the Chan- 
drasekhar limit. In fact, most WDs have masses peaked 
around O.6M , and most collisions are expected to have 

a total mass of 1.2 M Q , well below that limit. While 

IKatz & Pond d2012f) showed that the rate of such collisions 
may be as high as the rate of SNe la, it is unclear whether 
such collisions actually lead to explosions having the same 
observational properties as SNe la. 

Calculations of the thermonuclear expl osion of colliding 
WDs was attemp t ed by several groups dBenz et al.l 11989b 
Raskin et al. 2009; R osswog et al.ll2009t lLoren-Aguilar et al.l 
20 lOt iRaskin et al.l 120 10b IHawlev et al.H2012l) . using 3D hv- 
drodynamic simulations, motivated by the prospects of de- 
tecting rare SNe la in the cores of dense globular clusters. 
While the amount of 56 Ni synthesized in most of these sim- 
ulations was non negligible, the results were contradictory, 
yielding different amounts of 56 M. Furthermore, the ignition 
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site (contact region or shock region) is different between the 
different calculations for the same initial conditions. 

We use numerical simulations to calculate, from first 
principles, the explosion that is obtained from zero- 
impact parameter collisions of CO WDs with masses 
O.5,O.6,O.7,O.8,O.9M , covering the range of WD masses, 
including various equal and non-equal mass combinations. 
This problem is axis-symmetric, allowing the use of 2D nu- 
merical simulations, with high resolution (~Tew km cell size) 
that greatly supersede pre vious Eulerian appro aches (limited 
to > 100 km cell sizes) by IHawlev eTail (120121) . The cell vol- 
umes of our 2D calculations are smaller than the volumes of 
the particles in th e SPH runs (particl e mass divided by the 
highest density) by Raskin et al. (2010). Furthermore, we em- 
ploy ID planar toy models (see Figures [2] and |3]l to ensure 
that the ignition process, in the 2D simulation (see Figure [TJ 
is correctly resolved based on the fact that the curvature of the 
shock has little effect on the ignition process, which is largely 
a ID problem. Moreover, we used two differe nt 2D codes, 
Euleri an adaptive mesh refinement FLASH4.0 dDubey et al.l 
120091) and Arbitrary Eul erian Lagrang ian (ALE) code VUL- 
CAN2D (for details, see lLivndfl993l) which yield the same 
56 A// mass and the same location of detonation ignition. 

Numerical unstable burning occurs if the energy in a cell is 
significantly increased in a time shorter than the sound cross- 
ing time, t s = Ax/c s , where Ax is the length scale of the cell 
and c. 5 is the speed of sound. This is avoided by limiting the 
energy injection rate from burning, Q, to 0.1e/f„ where e is 
the internal energy of the numerical cell. This is implemented 
by appropriate renormalization of all burning rates within a 
cell if the limit is superseded. As the numerical resolution in- 
creases the renormalization becomes less severe, guarantee- 
ing correct convergence while avoiding premature ignition. 
We emphasis that the limiter does not modify the induction 
time, which is set at lower temperatures, where the stabil- 
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ity criterion is automatically satisfied. Therefore, the limiter 
does not modify the ignition process. Correct convergence is 
conformed by the calculations preformed with the ALE code, 
VULCAN2D, which did not require the use of a limiter at suf- 
ficiently high resolution (see upper panel of Figure |3)- Previ- 
ous Eulerian calculations, which violated this numerical sta- 
bility criterion, erroneously ignited too early a t the contact 
region for several scenarios (lHawlev et al.ll2012l) . where they 
should have ignited later at the shock region. This is the main 
reason for the reported contradiction between the SPH and the 
Eulerian codes. In our calculations, collisions involving low 
mass WDs (< 0.7) ignite at the shock region, while higher 
mass WDs ignite at the contact region, in agreement with the 
SPH simulations dRosswog et alj|2009h iRaskin et alj|2010l) . 

The amount of 56 Ni, synthesized in the various collisions is 
provided in Table Q] and shown in Figure [4] The collisions of 
typical CO WDs in the mass range 0.5- IM® produce 56 Ni 
masses in the range 0.1 - 1M Q in agreement with inferred 
range of SNe la, and thus provides a natural explanation for 
these SNe. This is the main result of this paper. 

A comprehensive comparison of the observational features 
of SNe to the collision model requires a detailed modeling of 
the complicated radiation transfer of optical light and is be- 
yond the scope of this paper. Next, we show that the model 
succeeds in satisfying two non-trivial, robust observational 
tests, which are independent of the optical radiation transfer 
details, and allow a direct comparison of the model to the ob- 
servations. 

The late time (> 50 days after peak) bolometric luminosity 
equals the total instantaneous 7-ray energy deposition rate, 
which is calculated exactly from first principles, by numeri- 
cally solving the transport of 7-rays, which is dominated by 
Compton scatterings. We find that the escaping fraction of 
7-rays as a function of time is the same in all of our mod- 
els, and agrees with each of the 24, low extinction (Galac- 
tic+host E(B-V)<0.3), Lboi/Mse^,- ob served in the sample of 
bolometric light curves compiled by iStritzinger et al.l (12005b 
(see Figures [S]|6]l. For illustration, the light curve of a simple 
Chandrasekhar toy model, with an exponential distribution of 
density (p oc e" v ) is shown. As can be seen, for 56 M mass 
of 0.15M Q (observed in some SNe), the Chandrasekhar light 
curve is higher than that from our model and also higher than 
the observations. This is due to the high optical depth of the 
massive 1.25M Q ejecta surrounding the 56 M. This ejecta is 
moving slightly slower (by ~ 10%) than that of the higher, 
0.8M©, 56 Ni models due to the lower amount of nuclear en- 
ergy released, and this leads to an additional ~ 20% increase 
in 7-ray deposition. 

The distribution of widths of the ( 56 M mass weighted) line 
of sight velocities for the collision models is compared to the 
widths inferred from nebular phase observations of SNe la 
in Figure [7] The model widths (v mo d) are obtained by fitting 
the LOS velocity distributions using a quadratic distribution 
(dMs6 Ni /dN oc l _ v 2 /v^ od ) commonly used and equivalent 
to homogenous expansion. The observational widths are ob- 
tained by fitting the spectra in the range 4800 -5700 A using 
the narrow spectra of 1991bg and 1999by as templates, by 
assuming local emission with a quadratic distribution of ve- 
locities. The amount of 56 M in the SNe is obtained by fitting 
intermediate bolometric light curves (f ~ 60 day) to the uni- 
versal injection function presented in Figure |5] which is well 

described by Ld ep0 sit = (1 +fe/40) 3 )~ 2 / 3 Ldecay, where td is the 
time from explosion in days and Ldecay is the energy released 



in 7-rays by Ni and Co. As can be seen, there is a clear 
correlation between the observed 56 Ni a nd the nebular phase 
velocity widths. This Mazzali relation (Mazzali et alJI 1998b 
is well reproduced by the collision model. The amount of 
SNe with nebular spectra and well described bolometric light 
curves is limited. In the top panel larger samples of SNe are 
used to show the continuous correlation of observational fea- 
tures of SNe la with 56 M yields. The strong correlation be- 
tween 56 Ni yields and WD mass in the collision model, Fig- 
ure |U is suggestive as the source fo r these correla tions, in- 
cluding possibly the Philips relation dPhillipsll 19931) . In fact, 
iPhillipsI d 1 993h had suggested a varying progenitor mass as a 
natural source for the relation. 

We emphasis that there are likely numerous detailed com- 
parisons of our calculations to observations. One very ro- 
bust and detailed prediction is the time dependent gamma-ray 
spectrum which can be calculated in a straight forward man- 
ner given the precisely known process of Compton scattering. 
While such observations may be very powerful, it is unclear 
when such observations will be available. Other directions 
which are closely related to the nebular velocity widths re- 
ported here, include the distribution of non-zero nebular line 
shifts, which are expected in the collision of non-equal mass 
WDs, as well as the unique line shapes resulting from the non- 
trivial velocity distribution in many of the collision scenarios. 

SNe la are used as standard candles, under the assumption 
that they would have no systematic change due to chemical 
evolution of the galaxy and physical evolution of the star- 
formation process, without any firm theoretical basis. Now 
for the first time we have a precise theory that allows us to 
calculate systematic changes if they exist, and the differences 
between ancient and modern SNe can be reliably estimated. 
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Table 1 

Summary of the results of the numerical simulations with FLASH, for the high resolution runs. Columns 1,2: The masses of the WDs (M l is the more massive 
WD); Column 3: the size of cells in the highest refinement level (we use 7 levels of refinement for these runs); Column 4: The type of ignition obtained in the 
simulation: for equal mass collisions the ignition is either at the shock region (s) or at the contact region (c). For non-equal mass collisions, the less massive star 
always ignites first, either near the shock region (s2) or near the contact region (c2), followed either by an ignition of the more massive star at the shock region 
(si) or no second ignition (xl); Column 5: The mass of 56 Ni synthesized in the explosion; Column 6: The total energy of the ejecta. 





*'*Z L J 


A.v [km] 


Ignition type 


56 Ni mass [M Q ] 


Btot [10 51 erg] 


0.5 


0.5 


9.5 


s 


0.10 


0.94 


0.55 


0.55 


9.0 


s 


0.19 


1.17 


0.6 


0.5 


9.5 


s2,sl 


0.19 


1.17 


0.6 


0.6 


8.5 


s 


0.29 


1.32 


0.64 


0.64 


8.1 


s 


0.37 


1.46 


0.7 


0.5 


9.5 


c2,sl 


0.27 


1.19 


0.7 


0.6 


8.5 


s2,sl 


0.44 


1.48 


0.7 


0.7 


7.6 


s 


0.50 


1.66 


0.8 


0.5 


9.5 


c2,xl 


0.28 


1.14 


0.8 


0.6 


8.5 


c2,sl 


0.36 


1.41 


0.8 


0.7 


7.6 


c2,sl 


0.52 


1.67 


0.8 


0.8 


6.8 


c 


0.73 


1.93 


0.9 


0.5 


9.5 


c2,sl 


0.64 


1.49 


0.9 


0.6 


8.5 


c2,xl 


0.52 


1.53 


0.9 


0.7 


7.6 


c2,sl 


0.48 


1.70 


0.9 


0.8 


6.8 


c2,sl 


0.72 


1.96 


0.9 


0.9 


6.1 


c 


0.76 


2.12 


1.0 


0.8 


6.8 


c2,xl 


0.89 


2.04 


1.0 


1.0 


5.3 


c 


1.25 


2.45 



density [g cm" 3 ] temperature [10 9 K] 

10 1 2 1 3 1 0" 10 5 10 6 1 7 0.7 1.7 2.9 7.0 




0.5 1 0.5 1 2 4 2 4 6 

r[10 lJ cm] r[10 9 cm] r [10 s cm] r [10 8 cm] 

Figure 1. The shock ignition process is illustrated using snapshots from the numerical calculation of a 0.64-0. 64 Mq collision done with FLASH with numerical 
resolution of ~ 4km. In panel (a) a logarithmic density map of the initial conditions of the calculation is shown. The black arrows indicate the direction of the 
velocity (assumed uniform) of each star. Following contact, a shock wave propagates from the contact surface toward the center of each star. The shock 
accelerates, leading to an increasing post shock temperature and r apidly increasing burning rates. Ignition of the detonation occurs once the induction time is 
shorter than the timescale for significant increase in burning rate (Zeldovich 1980). Panel (b) shows a density map at the time of ignition at the shock front 
(/ = 2.47 s). Two detonation waves begin to propagate from the ignition point, one inward and one outward. This process is fully resolved in our calculations, 
and the two ignition sites (one in each star) are clearly seen. The propagating detonation waves are evident in panel (d), which shows the same temperature 
map at a slightly later time, t = 2.57 s. Note that the mass surface where the ignition occurs becomes a contact discontinuity which is clearly seen in panel (d). 
An interesting feature seen in panels (c)-(d) is the outgoing jet, originating at the contact surface. This (planar) jet arrises due to the high pressures obtained at 
contact, which are not counter balanced from the outside. 
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Figure 2. The ignition process shown in Figure[T]is confirmed using a planar ID toy model evolved with Lagrangian and Eulerian schemes. This is possible 
since the hydrodynamical flow near the symmetry axis of the problem, where the ignition occurs, has approximate planar symmetry. The model consists of two 
equal mass WDs, with an initial planar density equal to the density on the axis of symmetry in the 2D model at the same distance from the contact surface. The 
gravitational field is mimicked by an adjustable constant acceleration (in time and space). In order to avoid artificial rapid early expansion due to the non physical 
tidal field, each mass element is forced by hand to fall at the constant external acceleration, until the point in time at which it is first shocked. The converged 
location of the ignition .Y„; n , obtained with a high resolution Lagrangian scheme is shown as a function of the adjustable acceleration, go: filled blue dots (shock 
region ignition) and red dots (contact region ignition). Beyond a critical acceleration of about 1.6 X 10 8 cms~ 2 the detonation occurs in the contact region. The 
acceleration go « 1.1 X 10 s , shown as a dashed line, approximately reproduces the 2D velocity profiles in Figure[T] The burning in the Eulerian code with 
resolution Aa ~ fewkm is unstable (see text), and leads to a premature detonation at the contact region (green cross) unless the burning is limited to be slower 
than the cell crossing time. The Eulerian scheme reproduces the correct ignition location when the limiter is included (orange cross). 
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Figure 3. The convergence of the 56 Ni mass as a function of resolution in our Eulerian 2D FLASH model is shown in black circles for the 0.64-0.64Mq 
collision simulation. The ALE V ULCAN2D code (red circle) reproduces the same value when c ompared at the same resolution. Our results are consistent with 
the SPH calculation performed by Raskin et al. 1 2010). In the Eulerian calculation performed by Hawlev et al. (2012), premature ignition at the contact surface 
occurs (blue point, see also Figure[2), due to burning which is faster than the cell sound crossing time. In the lower panel we show the decreasing error in the 56 Ni 
yield as a function of resolution in the ID Eulerian model. Based on the ID run, we estimate that the 5b Ni yield in our highest resolution FLASH ran (maximum 
level of refinement 8), is converged to about 10%. All other scenarios were calculated with lower maximum level of refinement 7, for which the 56 Ni masses 
are converged to about 20 — 30% accuracy (higher mass WDs have smaller radii and convergence is more easily obtained due to the smaller cell sizes for a given 
level of refinement). 
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Figure 4. The amount of 56 Ni synthesized in the numerical collisions is shown as a function of the mean mass of the colliding WDs at two resolutions and for 
equal and non-equal mass collisions. As can be seen, the yield function is converged and collisions of CO white dwarfs lead to the synthesis of ~ 0.1 - IMq of 
56 Ni covering the range of yields observed in the vast majority of SNe la. This is the main result of this paper. 
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Figure 5. Normalized bolometric light curves (black solid lines) of 24 SNe la, from the sample compiled by Stritzineer et al. (2005) . Panel a: Bolometric 
luminosity normalized to the peak luminosity. Panel b: Bolometric luminosity, normalized by the (time weighted) integrated luminosity, f dttL(t), which is 
proportional to the 5< Wi yield ( Katz et al. 2013). The integral is preformed to 80 days after the explosion, and the time of explosion is assumed to be 19 days 
prior to maximum. The instantaneous energy deposition by 7-rays in the expanding ejecta for all the collision models are shown (red dots), calculated using an 
exact Monte Carlo code, normalized by their time-weighted integral. At late times > 40day, the diffusion of optical radiation is short and this injected bolometric 
luminosity equals the emitted bolometric luminosity. The late time luminosity of the models is in excellent agreement with the observed light curves (see also 
a close up view in Figure|6j. The late time luminosity calculated from a simple Chandrasekhar model with 56 Ni yields of 0.8Mq (cyan dashed) and 0.15Mq 
(green dashed) is shown for comparison to illustrate that this agreement is far from trivial. Note that a much larger scatter is obtained when the light curves are 
normalized to the peak luminosity which is not directly proportional to the 56 Ni yield, as commonly believed. 
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Figure 6. Close up view of the bolometric light curves from Figure [5] 
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Figure 7. Bottom panel: Widths of the ( Ni weighted) line of sight (LOS) velocity distribution from the collision models (empty circles) compared to the 
widths inferred from nebular phase observations of SNe la (red points). Late time (>150 days) spectra are obtained from the Berkeley Supernova la Program 
(BSNIP) (Silverman et al. 2013), the Center for Astrophysics Supernova Program (Blondin et al. 2012), and the compilation from various sources by the Online 
Supernova Spectrum Archive (SUSPECT, http://suspect.nhn.ou.edu/~suspect/). The model widths (v rao[ j) are obtained by fitting the LOS velocity distributions 
using a quadratic distribution (dMs6 N J dN oc 1 — v z /v^ od ) commonly used and equivalent to an homologues expansion. The observational widths are obtained 
by fitting the spectra in the range 4800 — 5700 A using the narrow spectra of 1991bg and 1999by as templates, by assuming local emission with a quadratic 
distribution of velocities . The amount of 56 Ni in the SNe is obtained by fitting intermediate bolometric light curves (tj ~ 60day) to the universal injection 
function presented in Figure |5l which is well described by ^deposit = 

(i+fo/40) 3 r 2 / 3 z. 

decay where t(i is the time from explosion in days and Ldecay i s the energy 
released in 7-ray s by 56 Ni and 56 Co. As can be seen, there is a clear correlation between the observed s6 Ni and the nebular phase velocity widths. This so-called 
Mazzali relation i Mazzali et al. 1998), is well reproduced by the collision model. The number of SNe with nebular spectra and well described bolometric light 
curves is limited. In the top panel larger samples of SNe are used to show the continuous correlation of observational features of SNe la with 56 Ni yields. The 
strong correlatio n between 56 N i yields and WD mass in the collision model, Figure [4] is suggestive as the source for these correlations, including possibly the 
Philips relation (Phillips 1993). In fact, Phillips 1 1993) had suggested a varying progenitor mass as a natural source for the relation. 



